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Summary. We present a meshfree numerical solver for scalar conservation laws 
in one space dimension. Points representing the solution are moved according to 
their characteristic velocities. Particle interaction is resolved by purely local particle 
management. Since no global remeshing is required, shocks stay sharp and propagate 
at the correct speed, while rarefaction waves are created where appropriate. The 
method is TVD, entropy decreasing, exactly conservative, and has no numerical 
dissipation. Difficulties involving transonic points do not occur, however inflection 
points of the flux function pose a slight challenge, which can be overcome by a 
special treatment. Away from shocks the method is second order accurate, while 
shocks are resolved with first order accuracy. A postprocessing step can recover the 
second order accuracy. The method is compared to CLAWPACK in test cases and 
is found to yield an increase in accuracy for comparable resolutions. 
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1 Introduction 

Lagrangian particle methods approximate the solution of flow equations using 
a cloud of points which move with the flow. Examples are vortex methods [1], 
smoothed particle hydrodynamics (SPH) [13, 5], or generalized SPH methods 
[3]. The latter are typically based on generalized meshfree finite difference 
schemes [11]. An example is the finite pointset method (FPM) [9]. Moving 
the computational nodes with the flow velocity v allows the discretization of 
the governing equations in their more natural Lagrangian frame of reference. 
The material derivative = ■§[ + v ■ V becomes a simple time derivative. 
For a conservation law, the natural velocity is the characteristic velocity. In 
a frame of reference which is moving with this velocity, the equation states 
that the function value remains constant. Of course, this is only valid where 
the solution is smooth. In this case, characteristic particle methods are very 
simple solution methods for conservation laws. 
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In spite of the obvious advantages of particle methods, almost all numerical 
methods for conservation laws operate on a fixed Eulerian grid, even though 
significant work has to be invested to solve even a simple advection problem 
preserving sharp features and without creating oscillations. Leaving aspects of 
implementation complexity aside, two main reasons favor fixed grid methods: 
First, a fixed grid allows an easy generalization to higher space dimensions 
using dimensional splitting. Second, in particle methods one has to deal with 
the interaction of characteristics. While the former point remains admittedly 
present, the latter aspect is addressed in this contribution. 

Most methods which use the characteristic nature of the conservation law 
circumvent the problem of crossing characteristics by rcmeshing. Before any 
particles can interact, the numerical solution is interpolated onto "nicely" 
distributed particles, for instance onto an equidistant grid - in which case 
the approach is essentially a fixed grid method. The CIR-method [2] is an 
example. Such approaches incur multiple drawbacks: First, the shortest in- 
teraction time determines the global time step. Second, the error due to the 
global interpolation may yield numerical dissipation and dispersion. Finally, 
such schemes are not conservative when shocks are present. In practice, finite 
volume approaches, such as Godunov methods with appropriate limiters [14], 
or ENO [6J/WEN0 [12] schemes arc used to compute weak entropy solutions 
that show neither too much oscillations nor too much numerical dissipation. 

With moving particles, two fundamental problems arise: On the one hand, 
neighboring particles may depart from each other, resulting in poorly resolved 
regions. On the other hand, a particle may (if left unchecked) overtake a 
neighbor, which results in a "breaking wave" solution. The former problem can 
be remedied by inserting particles. The latter has to be resolved by merging 
particles. When characteristic particles interact (i.e. one overtakes the other) 
one is dealing with a shock, thus particles must be merged. 

In this contribution, we present a local and conservative particle manage- 
ment (inserting and merging particles) that yields no numerical dissipation 
(where solutions are smooth) and correct shock speeds (where they are not). 
The particle management is based on exact conservation properties between 
neighboring particles, which are derived in Sect. 2. In Sect. 3, we outline our 
numerical method. The heart of our method, the particle management, is de- 
rived in Sect. 4. There, we also show that the method is TVD. In Sect. 5, 
we prove that the numerical solutions satisfy the Kruzkov entropy condition, 
thus showing that the solutions we find are entropy solutions for any convex 
entropy function. In Sect. 6, we apply the method to examples and compare 
it to traditional finite volume methods using CLAWPACK. In Sect. 7 we 
present how non-convex flux functions can be treated. Finally, in Sect. 8, we 
outline possible extensions and conclusions. These include applications of the 
ID solver itself as well as possible extensions beyond the ID scalar case. 
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2 Scalar Conservation Laws 

Consider a one-dimensional scalar conservation law 

+ (/(w))^ = 0, u{x, 0) = uo{x) (1) 

with /' continuous. As long as the solution is smooth, it can be obtained by 
the method of characteristics [4]. The function u{x{t), t) is constant along the 
characteristic curve 

=x(0)+/'Kx(0),0))t . (2) 

For nonlinear functions / the characteristic curves can "collide" , resulting in 
a shock, whose speed is given by the Rankine-Hugoniot condition [4]. Dis- 
continuities are shocks only if the characteristic curves run into them. Other 
discontinuities become rarefaction waves, i.e. continuous functions which at- 
tain every value between the left and the right states. If the flux function / is 
convex or concave between the left and right state of a discontinuity, then the 
solution is either a shock or a rarefaction. If /" switches sign between between 
the two states, then a combination of a shock and a rarefaction occur. These 
physical solutions are defined by a weak formulation of (1) accompanied by 
an entropy condition. 

2.1 Conservation Properties 

Conservations laws conserve the total area under the solution 

oo 

d f 

— / u{x,t)dx = . (3) 

— oo 

The change of area between two moving points bi{t) and 62(f) is given by 

rb2{t) Mt) 



d_ 
dt 



u{x,t)dx = b'2it)u{b2{t),t) -b[{t)u{bi{t),t) + / ^((a;,*) da; 

6i(t) Jbi{t) 
= ib'2{t)u{b2{t),t)-b[{t)u{biit),t)) - if{u{b2{t),t))-fiuibi{t),t))) . 



If xi{t) and X2{t) are characteristic points, that is, points following the char- 
acteristics of a smooth solution as in equation (2), we have that x[{t) = f'{ui) 
and x'2{t) = f'{u2)- Therefore, the change of area between xi and X2 is 

{nU2)U2 - f{u,)u,) - {f{U2) - f{u,)) = [nu)u - fiu)]^ , (4) 

where [^(a;)]^ = g{b) — g{a). Equation (4) implies that the change of area 
between two characteristic points does not depend on the positions of the 
points, only on the left state ui and right state U2 and the flux function. Since 
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the two states do not change as the points move, the area between the two 
points changes Unearly, as does the distance between them: 

^^ix2{t) - x,it)) = - x[it) = f'M - r{m) = [f'iuZl . (5) 

If the two points move at different speeds, then there is a time to (which may 
be larger or smaher than t) at which they have the same position. Thus at 
time t = to the distance between them, and the area between them equal zero. 
From (4) and (5) we have that 

/ u{x,t)dx={t-to)-[nu)u-f{u)Zl , 

x2{t)-x,{t) = {t-to)-[nu)]:i . 

In short, the area between two Lagrangian points can be written as 

,-X2(t) 

u{x,t)dx = {x2{t) - xi{t))a{ui,U2) , (6) 



/ 



lxi(t) 

where a{u\,U'2) is the nonlinear average function 

[nu)u-f{u)]z i::f"{u)udu 

The integral form shows that a is indeed an average of u, weighted by /". 
This last observation needs one additional assumption: that the points xi and 
X2 remain characteristic point between t and Iq. That is, that a shock does 
not develop between the two points before to. Our numerical method relies 
heavily on the nonhnear average o(-, •)• 

Lemma 1. Let f he strictly convex or concave in [ul,Uu], that is f" < or 
/" > m {ul, uu). Then for all u\, U2 £ [ul, ujj], the average (7) is. . . 

1. the same for f and —f. Thus we assume f" > WLOG; 

2. symmetric, a{ui,U2) — a(u2,ui). Thus we assume Ui < U2 WLOG; 

3. an average, i.e. a{ui,U2) € {ui,U2), for ui ^ U2; 

4. strictly increasing in both ui and U2; and 

5. continuous at ui = U2, with a{ui,ui) = ui. 

Proof. We only include here the proof of 4. We show that a(ui, U2) is strictly 
increasing in the second argument. Let m < U2 < Us, Ui € [ul,Uu]. Then 

, , i::nu)udu+j:^\f"iu)udu 

«(^-"^) = JZrwu 

^ i::f"{u)du = "^^^'^^^ • 

Similar arguments show the result for the first argument. □ 
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3 Description of the Particle Method 

The first step is to approximate the initial function uq by a finite number of 
points xi, . . . , Xm with function values ui, . . . , Um- A straightforward strat- 
egy is to place xi, . . . , x„i equidistant on the interval of interest and assign 
Ui = uo{xi). More efficient adaptive sampling strategies can be used, since 
our method does not impose any requirements on the point distribution. For 
instance, one can choose Xi and Ui to minimize the error, using the specific 
interpolation introduced in Sect. 4. This strategy is the topic of future work. 
The points are ordered so that Xi < ■ ■ ■ < Xm ■ The evolution of the solution 
is found by moving each point Xi with speed f'{ui). This is possible as long 
as there are no "collisions" between points. Two neighboring points Xi{t) and 
Xi+i{t) collide at time t + Ati, where 

A positive Ati indicates that the two points will eventually collide. Thus, 
t + Ats is the time of the next particle collision^ , where 

Ats = min{Ati\Ati > 0} . 

i 

For any time increment At < Ats the points can be moved directly to their 
new positions Xi + f'{ui)At. Thus, we can step forward an amount Atg, and 
move all points accordingly. Then, at least one particle will share its position 
with another. To proceed further, we merge each such pair of particles. If the 
collision time Ati is negative, the points depart from each other. Although 
at each of the points the correct function value is preserved, after some time 
their distance may be unsatisfyingly large, as the amount of error introduced 
during a merge grows with the size of the gaps in the neighboring particles. To 
avoid this, we insert new points into large gaps between points before merging 
particles. In Sect. 4.1 we derive positions and values of the new particles that 
assure that the method is conservative, TVD, and entropy diminishing. 



4 Interpolation and Pcirticle Management 

The movement of the particles is given by a fundamental property of the 
conservation law (1): its characteristic equation (2). We derive particle man- 
agement to satisfy another fundamental property: the conservation of area 
(3). Using the conservation principles derived in Sect. 2, the function value 
of an inserted or merged particle is chosen, such that area is conserved ex- 
actly. A simple condition on the particles guarantees that the entropy does not 
increase. In addition, we define an interpolating function between two neigh- 
boring particles, so that the change of area satisfies relation (4). Furthermore, 
this interpolation is an analytical solution to the conservation law. 

^ If the set {i|Afi > 0} is empty, then Atg = oo. 
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4.1 Conservative Particle Management 

Consider four neighboring particles located at xi < X2 < X3 < X4 with associ- 
ated function values ui, U2, W3, U4. Assume that the flux / is strictly convex 
or concave on the range of function values [mini(ui), maxi(?ii)]. If U2 7^ M3, 
the particles' velocities must differ f'{u2) 7^ f'ius), which gives rise to two 
possible cases that require particle management: 

• Inserting: The two particles deviate, i.e. /'(U2) < /'(^a). If the distance 

X3 — X2 is larger than a predefined maximum distance d^ax^ we insert a 
new particle {X23, U23) with X2 < X23 < X3 and U23 chosen so that the area 
between X2 and X3 is preserved by the insertion: 

(0:23 - X2) a{u2, U23) + {x3 - X23) a(u23, U3) = {x3 - X2) 0(^2, M3) . (9) 

This condition defines a function, connecting {x2,U2) with {x3,U3), on 
which the new particle has to lie. This function is the interpolation defined 
in Sect. 4.2 and illustrated in Fig. 2. 

• Merging: The two particles collide, i.e. /'(W2) > /'(^s)- If distance 
X3—X2 is smaller than a preset value dmin (<^mm = is possible), we replace 
both with a new particle (X23, M23). The position of the new particle X23 is 
chosen with X2 < X23 < X3 and U23 is chosen so that the total area between 
Xi and X4 is unchanged: 

{X23 - Xi) a{ui,U23) + ix4 " a;23) a(M23: U4) 

= {X2-Xi)a{ui,U2) + iX3-X2) a,{u2,U3) + {X4-X3) a{u3,U4) . (10) 

Any particle (x23,W23) with X2 < X23 < X3 that satisfies (10) would be a 
valid choice. We choose 0:23 = ^^±23., and then obtain U23 such that (10) 
is satisfied. Figure 1 illustrates the merging step. 

Observe that inserting and merging are similar in nature. Conditions (9) and 
(10) for U23 are nonlinear (unless / is quadratic, see Remark 1). For most cases 
^23 = 2^2±m jg good initial guess, and the correct value can be obtained by 
a few Newton iteration steps. The next few claims attest that we can find a 
unique value U23 that satisfies (9) and (10). 

Lemma 2. The function value U23 for the particle at X23 is unique. 

Proof. We show the case for merging. The argument for insertion is similar. 
Prom Lemma 1 we have that both a{ui, •) and a(-, U4) are strictly increasing. 
Thus, the LHS of (10) is strictly increasing, and cannot be the same for 
different values of ■U23- □ 

Lemma 3. If X2 = X3 = X23, there exists U23 € [^2,^3] which satisfies (10). 

Proof. WLOG we assume that U2 < U3. First, we define 
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A = {x2 - xi) a{ui,U2) + {x4 - X2) a{u3,U4), and 
B{u) = {x2 — xi) a(Mi, u) + {x4 — X2) a{u, U4) . 



Equation (10) is now simply B{u23) = A. The monotonicity of a implies that 



Since a is continuous, so is B, and the existence of W23 follows the intermediate 



Corollciry 1. If particles are merged only according to Lemma 3, then the to- 
tal variation of the solution is either the same as before the merge, or smaller. 

Merging points only when X2 = X3 can be too restrictive. Fortunately, the 
following claim allows for a little more freedom. 

Theorem 1. Consider four consecutive particles (a;i,Uj) Vi = 1, . . . ,4. Merg- 
ing particles 2 and 3 so that X2i = ^^±23. yields U23 € [■U2,'M3] if 



Here the min and max of /" are taken over the maximum range of Ui, . . . , ^4. 
This condition is naturally satisfied if X2 = X3. 

Proof (outline). The full proof will be given in a future paper. The idea is to 
merge in two steps: First, we find a value u such that setting U2 — = u 
(while leaving X2 and X3 unchanged) results in the same area. Then, we merge 
the two points to U23. In the first step we bound u away from U2 and M3 (but 
inside [?i2,?i3]), and in the second step we bound |m23 — u\ from above. □ 

Theorem 2. The particle method can run to arbitrary times. 

Proof. Let ul = min^ Ui, uu = max^ u^, and C = max^uL-uu] \.f"{^)\'{^u—UL)- 
For any two particles, one has [/'(ui-i-i) — f'{ui)\ < C. Define Aa;^ ~ x^+i — 
Xi. After each particle management, the next time increment (as defined in 
Sect. 3) is at least Atg > . If wc do not insert particles, then in each 

merge one particle is removed. Hence, a time evolution beyond any given time 
is possible, since the increments Ats will increase eventually. When a particle 
is inserted (whenever two points are at a distance more than dmax), the created 
distances are at least , preserving a lower bound on the following time 
increment. □ 



B{u2) <A< B{u3) . 



(11) 



value theorem. 



□ 




(12) 
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Figure 1. Merging two particles 



Atl/'(ti2) 




A4i/'(<t) 



Ati/'(ni) 



Xl x{u) X2 Xgl, 



Figure 2. Definition of the interpolation 



4.2 Conservative Interpolation 

The particle management does not require an interpolation between points. As 
it stands, it complements the characteristic movcimciiit to yicJd a hill particle 
method for the conservation law (1) that can run for arbitrarily long times. 
Yet, for plotting the solution and interpreting approximation properties, it is 
desirable to define an interpolation that is compatible with the conservation 
principles of the underlying partial differential equation. We define such an 
interpolation between each two neighboring points {xi,Ui) and {x2, U2)- 

In the case ui = U2, we define the interpolation to be a constant. In the 
following, we describe the case ui ^ U2- Assume that / is strictly convex or 
concave in [ui,u2]. Therefore /'(w-i) / f'{u2)- Hence, as derived in Sect. 2.1, 
the solution either came from a discontinuity (i.e. it is a rarefaction) or it will 
become a shock. The time Ati until this discontinuity happens is given by (8). 
At time t + Ati the points have the same position x\ = X2 = Xsh_, as shown in 
Fig. 2. At this time the interpolation must be a straight line connecting the 
two points, representing a discontinuity at Xgh- We require any point of the 
interpolating function {x,u) to move with its characteristic velocity f'{u) in 
the time between t and t + Ati . This defines the interpolation uniquely as 

x{u) =x,-h if{u) - f\u,)) = XI + _ . (13) 

f'{U2) - f'iui) 

Defining function of u is in fact an advantage, since at a discontinuity 

characteristics for all intermediate values u are defined. Thus, rarefaction fans 
arise naturally if f'{u\) < f'{u2)- Since /" has no inflection points between 
ui and ?i2, the inverse function u{x) exists. However, it is only required at 
a single point for particle management. For plotting purposes we can always 
plot x{u) instead. 

Lemma 4. The interpolation (13) is a solution of the conservation law (1). 
Proof. Using that Xi{t) = /'(uj) for « = 1, 2 one obtains 



Solving Conservation Laws by Particle Management 



9 



Thus every point on the interpolation u{x, t) satisfies the characteristic equa- 
tion (2). □ 

Corollary 2 (exact solution property). Consider characteristic particles 

with xi{t) < X2 {t) < ■ ■ ■ < Xn (t) for t G [^1,^2]- At any time consider the func- 
tion defined by the interpolation (13). This function is a classical (i.e. con- 
tinuous) solution to the conservation law (1). In particular, it satisfies the 
conservation properties given in Sect. 2.1. 

Theorem 3 (TVD). With the assumptions of Theorem 1, the particle method 

is total variation diminishing. 

Proof. Due to Corollary 2, the characteristic movement yields an exact solu- 
tion, thus the total variation is constant. Particle insertion simply refines the 
interpolation, thus preserves the total variation. Due to Theorem 1, merging 
of particles yields a new particle with a function value U23 between the values 
of the removed particles. Thus, the total variation is the same as before the 
merge or smaller. □ 

Remark 1. The method is particularly efficient for quadratic flux functions. In 
this case the interpolation (13) between two points is a straight line, since /' is 
linear. Furthermore, the average (7) is the arithmetic mean a(ui, U2) = "^"^"'^ , 
since /" is constant. Consequently, the function values for particle insertion 
and merging can be computed explicitly. 

Rem,ark 2. The method has some similarity to front tracking by Holden et 
al. [7], and some fundamental differences. In front tracking, one approximates 
the flux function by a piecewise linear, and the solution by a piecewise constant 
function. Shocks are moved according to the Rankine-Hugoniot condition. In 
comparison, our method uses the wave solutions. Hence, in front tracking 
everything is a shock; in the particle method, everjd;hing is a wave. 



5 Entropy 

We have argued in Sect. 4.2 that due to the constructed interpolation the 
particle method naturally distinguishes shocks from rarefaction fans. In this 
section, we show that the method in fact satisfies the entropy condition 

7i{u)t + q{u) < (14) 

if a technical assmnption on the resolution of shocks is satisfied. We consider 
the Kruzkov entropy pair ??/c(w) = \u — k\, qkiu) = sign(u — k){f{u) — f{k)). 



10 



Yossi Farjoun and Benjamin Seibold 



As argued by Holden et al. [8], if (14) is satisfied for r]k, then it is satisfied 
for any convex entropy function. Relation (14) implies that the total entropy 
/ r]k{u{x))dx does not increase in time for all values of k. Using the inter- 
polation (13) we show that the numerical solution obtained by the particle 
method satisfies this condition. 

Lemma 5 (entropy for merging). Consider four particles located at x\ < 
X2 = Xz < X4, with the middle two to be merged. We consider the case f" > 0, 
i.e. U2 > U3 WLOG.^ If the resulting value U23 satisfies ui > U23 > U4, then 
the Kruzkov entropy does not increase due to the merge. 

Proof. We consider the segment [a;i,a;4]. Let u{x) and u{x) denote the inter- 
polating function before resp. after the merge. The area under the function 
is preserved. We present the proof for k < M23. For k > U23 the proof is sim- 
ilar. The interpolating function u is monotone in the value of its endpoints, 
thus u{x) < u{x) for x E [x2,x4]. Since = a; — 20(— a;), where Q{x) is the 
Heaviside step function, we can write 

/'*^4 P*^4 



The assumption of Lemma 5 implies that shocks must be reasonably well 
resolved before the points defining it are merged. It is satisfied if left and 
right of a shock points are not too far away. In the method, it can be ensured 
by an entropy fix: A merge is rejected a posteriori if the resolution condition 
is not satisfied. Then, points are inserted near the shock, and the merge is 
re-attempted. It remains to show in future work that with this procedure 
Theorem 2 still holds. 

Theorem 4. The presented particle method yields entropy solutions. 

Proof. During the characteristic movement of the points the entropy is con- 
stant, since due to Corollary 2 the interpolation is a classical solution to the 
conservation law. Particle insertion does not change the interpolation, thus 
it does not change the entropy. Merging does not increase the entropy if the 
conditions of Lemma 5 are satisfied. □ 



^ For the case /" < 0, all following inequality signs must be reversed. 




\u—k\ dx . 



Thus, the entropy does not increase due to the merge. 



□ 
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6 Numerical Results 

The particle method is particularly well suited for initial conditions that are 
composed of similarity solutions. By construction, the movement of the parti- 
cles yields the exact solution as long as the solution is smooth. General initial 
conditions can be approximated by the interpolation (13). Good strategies of 
sampling initial conditions shall be addressed in future work. Figure 3 shows 
a smooth initial function uo{x) = cxp {—x'^) cos (ttx) and its time evolution 
under the flux function f(u) = ju^. The curved shape of the interpolation 
is due to the nonlinearity in /'. At time t = 0.25 the solution (obtained by 
CLAWPACK using 80000 points) is still smooth, and thus represented exactly 
on the particles. At time t = 8 shocks and rarefactions have occurred and in- 
teracted. Although the numerical solution uses only a few points, it represents 
the true solution well. 




The accuracy of the particle method is measured numerically. We consider 
the flux function and initial conditions as used in Fig. 3. For a sequence of 
resolutions h, the initial data are sampled, and the particle method is applied 
{dmax = l-9/i). Figure 4 shows the L^-crror to the correct solution (obtained 
by a computation with much higher resolution, verified with CLAWPACK). 
While the solution is smooth {t = 0.25), the method is second order accu- 
rate, as is sampling the initial data. After a shock has occurred {t = 0.35), 
the approximate solution (dots) becomes only first order accurate, since the 
shock has just been treated by particle management, thus an error of the or- 
der height X width of the shock is made. A postprocessing step (squares) can 
recover the second order accuracy: At merged particles, discontinuities are 
placed so that the total area is preserved. 



7 Non- Convex Flux Functions 

So far we have only considered flux functions with no inflection points (i.e. /" 

has always the same sign) on the region of interest. In this section, we gen- 
eralize to flux functions for which /" has a finite number of zero crossings 
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10^ 10' 10"^ 10' 



resolution h resolution h 

Figure 4. Error to the correct solution before and after a shock 

< • • • < M^. Between two such points u € [u*, u*_^_-^] the flux function is ei- 
ther convex or concave. We impose the following requirement for any set of 
particles: Between any two neighboring particles for which /" has opposite 
sign, there must be an inflection particle {x, u*). Thus, between two neighbor- 
ing particles, / has never an inflection point, and most results from the pre- 
vious sections apply. The interpolation between any two particles is imiquely 
defined by (13). It has infinite slope at the inflection points, but this is harm- 
less. The characteristic movement of particles is the same as for flux func- 
tions without inflection points. The only complication is merging of particles 
when an inflection particle is involved: The standard approach, as presented 
in Sect. 4.1, removes two colliding points and replaces them with a point of a 
different hmction value. If an inflection particle is involved in a collision, we 
must merge points in a different way so that an inflection particle remains. 




We present one such special merge for dealing with a single inflection 
point (we do not consider here the interaction of two inflection points). Also, 
we only consider a collision where the positions arc exactly the same. Since the 
inflection particle must remain (although its position may change), we consider 
five neighboring particles and not four as before. Let {xi,Ui),i = 1, . . . , 5 be 
these particle so that X2 — X3, f"{u^) = 0, and (WLOG) /"' > 0, i.e. the 
inflection particle is the slowest. The other cases are simple symmetries of this 
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situation. We present three successive steps to finding the final configuration 
of the particles. Each next step is attempted if the previous one failed. 

1. Remove particle 2 and increase to satisfy the area condition. This fails 
if X3 needs to be increased beyond X4. 

2. Remove particle 2, set X3 = X4 and increase both to satisfy the area con- 
dition. This fails if X3,a;4 need to be increased beyond X5. 

3. Remove particle 4, set X3 = X5 and find U2 to satisfy the area condition. 
This cannot fail if the previous two have failed. 

Formally, the resulting configuration could require another, immediate, merge 
(since X3 — X4 or X3 = X5). However, we need not merge these points as they 
move away from each other. The five point particle management guarantees 
that in each merging step one particle is removed, thus Theorem 2 holds. 




1 2 3 1 2 3 1 2 3 



Figure 6. Numerical results for the Buckley-Leverett equation 



As numerical evidence of the performance, we apply our method to the 
Buckley-Leverett equation (see LeVeque [10]), defined by the flux function 
f{u) — 3 1" -3 . It is a simple model for two-phase fluid flow in a porous 
medium. We consider piecewise constant initial data with a large downward 
jump crossing the inflection point, and a small upward jump. The large jump 
develops a shock at the bottom and a rarefaction at the top, the small jump is 
a pure rarefaction. Around t = 0.2, the two similarity solutions interact, thus 
lowering the separation point between shock and rarefaction. Figure 6 shows 
the reference solution (solid line, by CLAWPACK using 80000 points). The 
solution obtained by our particle method (dots) is compared to a second order 
CLAWPACK solution (circles) of about the same resolution. While the finite 
volume scheme loses the downward jump very quickly, the particle method 
captures the behavior of the solution almost precisely. Only directly near the 
shock inaccuracies are visible, which are due to the crude resolution. The 
solution away from the shock is nearly unaffected by the error at the shock. 
Note that although we impose a special treatment only at the inflection point, 
the switching point between shock and rarefaction is identified correctly. 
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8 Conclusions and Outlook 

We have presented a particle method for ID scalar conservation laws, which 
is based on the characteristic equations. The associated interpolation yields 
an analytical solution wherever the solution is smooth. Particle management 
resolves the interaction of characteristics locally while conserving area. Thus, 
shocks are resolved without creating any numerical dissipation away from 
shocks. The method is TVD and entropy decreasing, and shows second-order 
accuracy. It deals well with non-convex flux functions, as the results for the 
Buckley-Leverett equation show. The particle method serves as a good alter- 
native to fixed grid methods whenever ID scalar conservation laws have to be 
solved with few degrees of freedom, but exact conservation and sharp shocks 
are desired. An application, which we plan to investigate in future work, is non- 
linear flow in networks (e.g. traffic flow in road networks). For large networks, 
only a few number of unknowns can be devoted to the numerical solution on 
each edge. We regard the current work as a first step towards a more general 
particle method. Future work will focus on three main directions: 

• Source terms: Source terms in an equation ut + {.f{u))^ = g{x,u) could 
be handled using a fractional step method: In each time step, we would 
first move the particles according to Ut + {f{u))^ = (including particle 
management), then change their values according to an integral formula- 
tion of Ut = g{x,u). In the latter step, the constructed interpolation can 
be used. 

• Systems of conservation laws: The particle method is based on similar- 
ity solutions of the conservation law. For simple systems, such as the shal- 
low water equations in ID, the analytical solutions to Riemann problems 
are known. Two complications arise in the generalization of the method: 

- To connect two general states in a hyperbolic system, intermediate 
states have to be included. 

- For a general system it is not clear at which velocity to move the 
particles. 

• Higher space dimensions: Scalar conservation laws in higher space di- 
mensions can be reduced to ID problems to be solved by fractional steps. In 
principle, this dimensional splitting can be used with the particle method. 
However, remeshing would be required between the different spacial di- 
rections, thus the benefits of the meshfrec approach would be lost. For 
the generalization to a fundamentally meshfree approach in higher space 
dimensions, the following problem has to be overcome: In ID one is never 
truly meshfrec, since the points have a natural ordering. The method uses 
this in the interpolation and to detect shocks. In 2D/3D shocks can occur 
without particles colliding, as they can move past each other. Other mesh- 
free methods, such as FPM applied to the Euler equations, circumvent this 
issue by using the pressure to regulate shocks. 
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